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Abstract 



Lattice NRQCD with leading finite lattice spacing errors removed is 
used to calculate decay constants of mesons made up of heavy quarks. 
Quenched simulations are done with a tadpole improved gauge action contain- 
ing plaquette and six-link rectangular terms. The tadpole factor is estimated 
using the Landau link. For each of the three values of the coupling constant 
considered, quarkonia are calculated for five masses spanning the range from 
charmonium through bottomonium, and one set of quark masses is tuned 
to the B c . "Perturbative" and nonperturbative meson masses are compared. 
One-loop perturbative matching of lattice NRQCD with continuum QCD for 
the heavy-heavy vector and axial vector currents is performed. The data are 
consistent with afy cx \JMy a and Jb c = 420(13) MeV. 

PACS number(s): 12.38.Gc, 12.39.Jh, 13.20.Gd, 13.20.Jf 
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I. INTRODUCTION 



The gross features of quarkonia are well described by quenched lattice NRQCD 
H|-0|. However, spin splittings between vector and pseudoscalar mesons tend to be 
underestimated — especially when relativistic corrections are included . Using the quark 
model, spin splittings are due to a short-ranged Fermi-Breit interaction and essentially mea- 
sure the square of the meson's wave function at the origin. Since they are not well understood 
within lattice NRQCD it is of interest to investigate physical quantities which probe the same 
physics. This is provided by the mesonic decay constants which, in the quark model, are 
proportional to the wave function at the origin. 

In this paper, the focus is on the vector decay constants of quarkonia and the pseudoscalar 
decay constant of the B c meson. The calculated quarkonia constants can of course be 
compared with experiment — for bottomonium this provides a test of the wave function at 
the origin, which is not possible from the spin splitting since the rjb has not been observed 
yet. On the other hand, the B c has only recently been discovered and its decay constant 
is not measured; however, comparison between our results and previous lattice NRQCD and 
other model predictions can be made. 

A previous study |§ of the vector decay constant found large corrections from the per- 
turbative matching between the continuum and lattice currents. Another study reported 
a rather imprecise value for simulations which included the order v 2 relativistic corrections 
of the vector current. Previous studies [TOP] of the B c decay constant reported results 



which did not include matching or relativistic corrections in the current itself. In this paper 
we remove the leading finite lattice spacing errors in the fermion and gauge actions in a 
symmetric fashion and perform (1) more precise simulations with the inclusion of the or- 
der v 2 relativistically corrected currents and (2) one-loop matching between the lattice and 
continuum currents at lowest order in v 2 . 

The paper is arranged as follows. In Sec. |IJ lattice NRQCD is introduced and tadpole 
improvement is discussed. Continuum and lattice matrix elements are matched up via 
perturbation theory in this paper. Thus, in Sec. [1 1 1] lattice perturbation theory is set up. 
Appendix [A| lists the free gluon propagator and the integrand of the amputated axial vector 
correction. A discussion of the matching with continuum QCD is presented in Sec. [TV 
(Appendix |B] shows the origin of the '1/V terms in the continuum QCD result). A general 
discussion of lattice decay constants, along with our conventions, is presented in Sec. [V[ and 



finally our results are discussed in Sec. VI 



II. LATTICE NRQCD 

The fermion NRQCD Lagrangian is discretized in a symmetric fashion with the leading 
finite lattice spacing errors in the spatial and temporal derivatives removed. All links are 
tadpole improved by dividing by it , the average link in the Landau gauge. The choice 
of 'Landau link' over 'average plaquette' tadpole is made based on spin splitting studies 
of quarkonia [||] where improved scaling behavior was observed. Tadpole improvement in 
general improves the matching between perturbation theory and lattice simulations [12[] (see 
Table | for our evidence for this) . The fermion Lagrangian we use is 
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t t/i aSH\ aHz\ n U\ t -\(, aH \" a5H s 



2n J t uo V 2n /t-i V 2 



t-i 



A( 2 ) X[7 _a 2 A( 4 ) a(A( 2 ))' 



#o = , SH = i '— , (lb) 

2m 24m Vonm 1 

u = (^ReTrU^ , 9^ = 0. (lc) 

n is the stability parameter chosen to satisfy n > 3/ (ma). A^ 2 ) is the gauge-covariant lattice 
Laplacian, and A^ 4 ) is the gauge-covariant lattice quartic operator (I]j-D 4 )|] 

The above form of the fermion Lagrangian is convenient for defining Feynman rules. In 
lattice simulations the action leads to an evolution equation for the quark propagator of the 
form 

aSH\ / aH \ n U\ t _ x ( aH \ n / a5H\ 



a M - ( 1 - {1 - ^ {1 - ^ L {1 - — G(x,t - 1) 

a5\x) , (2) 



where in this equation the source has been placed at the origin for convenience and 
G(x,t) = for t < 0. 

The gauge action is tadpole improved with leading finite lattice spacing errors removed 
by six-link rectangles ||13||: 



S G = f3 ^Re Tr (1 - U pl ) - £ ^ReTr (1 - U rt ) . (3) 

pi 6 ZUU rt 6 

In lattice NRQCD meson masses can not be calculated directly from meson time correla- 
tion functions. However, using meson correlators at zero and nonzero momentum, a 'kinetic' 
mass is obtained by 

E(P) ~ 25(0) = , (4) 

where E(p) is the simulation energy and the momentum p is expressed in units of 2ir/(Na), 
where 'TV a' is the spatial extent of the lattice. 



1 For the standard definitions of these lattice derivatives see for example Ref. 
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III. PERTURB ATIVE LATTICE NRQCD 



The Feynman rules are derived from the actions of Eqs. (|I|) and (|3|) by making the 
replacement 



exp^agAVx)^ 



(5) 



and expanding in g. This is not quite true because a gauge fixing term must be added to 
the gauge action before expanding in g. The standard covariant gauge-fixing term in lattice 
perturbation theory is 



(6) 



where the gluon field is to be expanded about the midpoint of the link: 

m d 4 q 



KM 



-A a Jq)exp[iq-(x + ti/2)\ . 



U(2tt) 4 ^ 

To be complete, the fermion field in momentum space is defined by 

^ d 4 p 



(7) 



ip(x) 



tt (2tt) 



-M>{p) exp (ip ■ x) 



(8) 



The gluon propagator follows from the quadratic part of Sq + Sgf- See Appendix [A] for its 
explicit form. Note that we use the Feynman gauge (£ = 1) for perturbative lattice NRQCD 
self-energy and vertex corrections in this paper. 

The coupling g in this perturbative identification [Eq. (|5|)] that we use is the so called 
"boosted coupling." It is defined by rewriting the gauge action as 



S, 



6 

g 2 



yiReTrfl-^ - -V-ReTrfl 



(9) 



The 5/3 and —1/12 are necessary to match the continuum action with its usual normaliza- 
tion. Comparing Eqs. (H) and (0) gives 
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a 



4tt 3 



(10) 



The benefits of using this boosted coupling can be seen by comparing the Landau link 
calculated nonperturbatively and perturbatively [0: Define u^ 1 by 



Wo 



a u 



(2) 

5 



(11) 



where u is calculated nonperturbatively and a(j3, u ) is given by Eq. (|10|) . Then com- 
pare with a perturbatively calculated Uq\ This comparison is shown in Table the close 
agreement justifies the use of a boosted coupling. 
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Given the Feynman rules, the perturbative matching factors for the decay constants 
follow. We discuss this matching in detail in the next section. A "perturbative" meson mass 
can also be defined. For example in the equal mass case 

M pert = 2(mZ m — E ) + E sim , (12) 

where Z m is the mass renormalization and Eq is the energy shift of a quark, and E S i m is the 
simulation energy extracted from the meson correlator. For details on how Z m and Eq are 



defined see Ref. 14 



IV. MATCHING WITH CONTINUUM QCD 

Lattice NRQCD is an effective field theory that is fundamentally different from contin- 
uum QCD — the ultraviolet divergences in the matrix elements of interest are not the same. 
After renormalization, these differences are finite (since the infrared divergences are the 
same), but nevertheless in order to obtain continuum QCD results, a matching step must 
be performed. In this section we will focus on the axial vector current matching since the 
vector current case follows in a similar manner. Details of the axial vector current matching 
will be given, and in the end, results for the vector current case will also be shown. 

In general for the annihilation decay of the B c meson, matching leads to 

(0\hol 5 c\cb) = Z match (0|6 7 o75c|c5) + 0{v 2 ) , (13) 

cont lat 

where v is the relative velocity of the b and c. In performing this matching in this paper we 
neglect the order v 2 terms and calculate Z match in one-loop perturbation theory. For this 
system a ~ .25 and v 2 ~ .2: the expansion parameters are somewhat less than unity. 

Although the matching can be done in one step,0 we will use two steps for pedagogical 
reasons. The two step matching procedure starts with the matching of lattice NRQCD 
(1NRQCD) with continuum NRQCD (cNRQCD). Then continuum NRQCD is matched with 
continuum QCD (cQCD). In the end a one-loop formula, relating the simulated lattice 
NRQCD matrix element with its continuum counterpart of full QCD, is obtained. 

We will use the following notation: the amputated vertex correction will be written as 
i g 2 5V\ including the wave function renormalization factors, the total matrix element will 
be written as '1 + g 2 (8V + SZb/2 + 5Z c /2)' times the tree level amplitude; and we will use 
different regulators as given below — but our renormalization scheme is always the on-mass- 
shell scheme. 



A. Matching lattice NRQCD with continuum NRQCD 

The current matrix element is ultraviolet finite in NRQCD. However infrared divergences 
do arise, and in this subsection a gluon mass is used as a regulator. The matrix elements 
are calculated in the limit 



2 We performed the matching in one step also (matching continuum QCD directly with lattice 
NRQCD) and obtained the same result as what follows from the discussion below. 
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1 > A/ m red > v 







(14) 



where A is the gluon mass, v is the relative velocity and m red is the reduced mass of the 
system, 



l/m red = l/m b + l/m c . 



(15) 



In continuum NRQCD it is convenient to work in the Coulomb gauge. The transverse 
gluons coupling to the quarks is suppressed by powers of v , so only the Coulomb exchange 
term needs to be considered. The wave function renormalization factor is 



9 2 SZ 



cNRQCD 



-iCf9 2 J 



d 4 k 



1 



(27r)4(k2 + A2)(fc +£-^) 2 ' 



(16) 



where Cf = 4/3. Performing the & integration, we see that the result is trivial, 
3Z\cNrqcd = 0- The amputated vertex correction is 



9 2 6V 



d 4 k 



cNrqcd -"j (27r) 4 (k 2 + A 2 )(A; 
Performing the ko integration, we are left with 

9 2 SV 



(17) 



2mt 



-9 r 2 f 1 _ 2m red9 2 

cNRQCD ~ Zm ^f9 J ^3 ^ + ^2 " ^ ■ 



(18) 



In summary, for continuum NRQCD with the scales proportioned as in Eq. (|14D through 
one loop we have 



(0|x^c|c6> 



cNRQCD 



vkc 



1 + 



2m red g 2 
3ttA 



(19) 



where ip c and Xb are non-relativistic c and b fields respectively. 

Lattice NRQCD is defined in Euclidean space. However, since the zeroth component 
of the axial vector current is identical in Euclidean and Minkowski space, one can match 
directly with the continuum NRQCD result .0 The detailed form of the integrands is not 
given here (see Appendix [X]). Instead, we will just write the form of the final result 



{0\Xb*Pc\cb) 



INRQCD 



1 + 9 2 



2m, 



red , r-T7 



3ttA 



+ 5V lat + 5Zf at /2 + 8Z*J2 



(20) 



where the bar on 5Vi a t signifies the separation of the linear infrared divergence. The integra- 



tion routine VEGAS [J5j is used to perform the integrations. Only infrared finite integrands 
are integrated with VEGAS. We subtract and add low three-momentum versions (leaving the 
energy variable alone) of the original integrand, where the subtraction is chosen to (1) make 



3 For the vector case there is just a trivial factor of i that must be kept track of: 7 



-*7, 
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the original integrand infrared finite and (2) leave a simpler infrared divergent integrand 
that can be integrated analytically (or else repeat the subtraction procedure on this new 
term also). In an equation 



d 4 kl(k) 



d 4 k 



I(k) - Isub(k) 6(c 2 - k 2 



d A kI sub (k)9(c 2 -k 2 ) 



(21) 



where the step function is chosen for convenience in the analytic integration step- 
arbitrary parameter less than 'tt' (for example '1' or '2'). In this way, the linear infrared 
divergence is analytically separated as we have written it in Eq. (f2H|). 

To conclude, divide Eq. (|19|) by Eq. (^0|). Through one loop with the scales proportioned 
as in Eq. (|TJ) 



(0\Xbi>c\cb) 



cNRQCD 



(0\XbA\cb) 



INRQCD 



1 - 9 2 (SV lat + 8ZIJ2 + 8Z h lat /2 



(22) 



where the remaining matching factor is infrared finite; however, a logarithmic divergence 
cancels between the remaining terms. 

B. Matching continuum NRQCD with continuum QCD 

First, our conventions for this subsection: dimensional regularization regulates the in- 
frared (IR) and ultraviolet (UV) divergences;^ we will work with finite but small relative 
velocity v ; 75 is chosen to anticommute with all other gamma matrices; and as already 
mentioned, our renormalization scheme is the on-mass-shell scheme. 

The necessary matching calculation for this subsection has been performed already by 
Braaten and Fleming [ l6| . Details are omitted here and can be found in Ref. [Rj. However, 
one intermediate step of the calculation is discussed in Appendix |B], that is, the origin of 
the '1/V terms in the results. 

The one-loop renormalization of the axial vector current matrix element consists of am- 
putated vertex and wave function renormalization corrections. In continuum NRQCD it is 
convenient to use the Coulomb gauge. Here the wave function renormalization correction 
vanishes [see Eq. ([16])], and the amputated vertex correction is UV finite. However there is 
an IR divergence. The full result in continuum NRQCD is 



<0|X6^c|c6) 



cNRQCD 




2m red v 

— 2 log h log 477 — 7 

ClR ^ 1 



(23a) 



In continuum QCD it is convenient to use the Feynman gauge. The matrix element is UV 
divergent, but these divergences cancel between the amputated vertex and wave function 



4 The usual D = 4 — 2e identification is made. 
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renormalization factors. There is a real IR divergence that also cancels, but the final result 
has an imaginary IR divergence. The full result in continuum QCD is 



(0|&7o75c|c6) 



cQCD 



1 + 



g 2 | if_ 3 m b - m c 



log 



m b 



6tt 2 I v 2 m& + m c m, 



in ( I 2m red v 

— 2 log h log 47r — 7 



(23b) 



For the respective one-loop integrals to converge em < and ejjv > 0. Also, the diagrams 
are evaluated slightly above threshold, 



q 2 = (nib + ^c) 2 + mbm c v 2 + 0(v 4 ) 



(24) 



thus the b and c can simultaneously go on-mass-shell giving rise to the imaginary contribu- 
tions. Dividing these continuum QCD and NRQCD results [Eq. (|23|) 1, through one loop and 
for small v, gives for the final result 



(0|&7o75c|c&) 



cQCD 



{0\Xbtpc\cb) 



cNRQCD 



1 + 



g 2 /3m(,-m c 



r > \ o i lo S 

D7r z \2 nib + vn c m, 



(25) 



where 'fb7o75M c = Vb£c + 0(v 2 y has been used. 



C. Matching lattice NRQCD with continuum QCD 

In both of the previous two subsections the ultraviolet and infrared divergences canceled 
in the final result, thus the two steps can be combined. Substituting Eq. fl2"2] ) into Eq. 
gives for our final result 



(0|&7o75c|c&) 



cQCD 



INRQCD 



1 + 



g 2 (3m b -m c m b " 

n 1 log 3 

bn z \2 nib + m c i^c 



g 2 (5V lat + 5ZIJ2 + 5Z h lat /2) 



(26) 



which as already mentioned is the same result we obtained by directly matching lattice 
NRQCD with continuum QCD in one step with the scales proportioned as in Eq. (|i~4]). 

As promised, the result for the vector case will also be given. The procedure is the same 
as for the axial vector, and the result is very similar — take the previous with equal quark 
masses and '—3 — > —4': 



(0\QiQ\QQ) 



cQCD 



(o\xq^q\QQ) 



INRQCD 



1 + 



9 

6tt 2 



-4) - g 2 (5V lat + 5Z? at 



We will use the following notation to report these matching contributions: 

(0\J\M)\ 

cQCD ~ ^match 

(0\J\M)\ 

INRQCD ■ 

Tables [TT] and [TIT] show our results. 



• (27) 



(28) 
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V. DECAY CONSTANTS 



A meson propagator is given by 

G(p,t) = /^exp[-ip ■ (x-x )]J(x,t)J t (x ,t )^ 

= (j2 e *P[-ip- (x-xo)]Tr [(T x G XX0 )(G XX0 T X0 ) j ^j , (29) 

where the trace is over spin and color, and G xxo is the quark propagator of Eq. @. J(x, t) = 
XxF x ip x is a non-relativistic current with T x = Q x u) x , where Q x interpolates the meson of 
interest^ and u x is a smearing operator chosen in a gauge-covariant fashion: 

u x = [l + eA^(x)r . (30) 

We set e = 1/12 and tune the smearing parameter n s to maximize the overlap with the state 
of interest. The range 7-30 for n s was found to be sufficient, with the P-wave requiring 
about twice as much smearing as the S-wave and the smearing parameter increasing for 
decreasing quark mass. 

The asymptotic form of this non-relativistic meson propagator is 

G (P^) t^w^o |(0|J(0)|p>| 2 exp[-£(p)(t-t )]. (31) 

In continuum Minkowski notation, this current matrix element for a pseudoscalar and vector 
at rest is related to the respective decay constant by 

i f R Mr 

(OlfnwIBc) = 1 j=j^ , (32a) 

where e is a polarization vector and for matching purposes a non-relativistic norm is assumed. 

Simulations are performed with order v 2 relativistic corrections included in the currents. 
The relativistically corrected interpolating operator for a non-relativistic vector meson is 
given by 

n " = " + ( A<2,v - + <^ • A<2, ) - ( CT ' A ') ffi (" ■ A ) ■ (33) 

where A is the gauge-covariant symmetric lattice derivative and A^ is the gauge-covariant 
lattice Laplacian. For the B c the relativistically corrected interpolating operator is 

A(2) 

Mb! = 1 + 5-5- • (34) 



5m 



red 



b Q x = I, <t, and A for the l S$, 3 5i, and l P\ states respectively. 
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Our final decay constants are reported using the following notations: Let 'src' represent 
L (local) or S (smeared). 'X' is l P\, PS (pseudoscalar), or V (vector). Three correlators are 
of interest: 

(1) C L X = ^4(x,t)4 t (x ,t )^ = Z x exp [-E x (t-t ] 

(2) C s x = ^4(x,t)j| t (x ,t )^ ee Z s x exp [-E s x (t-t ] 
and (3) C S / el = (^J^ l ( x ,t)4\x ,t )) = Z^ 



(35) 



where J x rel (x., t) = Xx^x^x and the right-hand-side of these equations assumes t ^> to- 
Given this, we report two forms of decay constants: 

fxM x 



v/2M 



o 2 Zmztchy %x 



X 



and 



f r x el M x 



V2M 



a 2 -^match" 



X 



I ryL yS,rel 
^X ^X 



(36) 



(37) 



where f x l includes the order v 2 relativistic corrections in the currents and Z ma tch is the 
perturbative matching factor of Eq. (ESI). The spirit of these definitions is from Ref. [17]. 



VI. RESULTS AND DISCUSSION 

Tables [TV] and [V] show the simulation parameters and compare M pert and M^ in . The 
lattice spacing is fixed by the splitting between S- and P-wave states which is taken to be 
458 MeV for all simulations. For each coupling, five quark masses spanning the range from 
charmonium through bottomonium are used, and one set of quark masses is tuned to the 
B c which we take to be 6.35 GeV for the spin averaged ground state. 

The data sample includes 1200 quenched gauge field configurations at /3 — 7.4 (10 3 x 
16), and 1600 configurations at (3 — 7.2 and 7.3 (8 3 x 14). A standard Cabbibo-Marinari 
pseudo heat bath is used to generate the gauge-field configurations. After 4000 thermalizing 
sweeps, the number of updates between measurements is 30. Autocorrelation times for 
the correlation functions are checked and satisfy r < |. Multiple sources (N/2) along the 
spatial diagonal are used to measure the local-smeared meson correlators with the number 
of smearing steps optimized to have the best overlap with the state of interest. Single 
exponential fits to the correlation functions are used to get the best estimates of the masses 
and decay constants. Effective mass plots are generated and used to choose the interval in 
which to do the fits. In all cases acceptable Q values are obtained. All statistical errors are 
estimated by the bootstrap method with twice as many bootstrap ensembles as there are 
configurations. As expected, simulation energies from L-L, L-S, and Lrel-S correlators are 
all consistent with each other. 

Fig. [I] shows our main result for quarkonia: 

af v oc ^M v a . (38) 
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For a tabular form of this information see Table [VJ| . This square-root dependence is quite 
interesting. Shortly after the discovery of charm, Yennie ]n| noticed this same dependence 
from empirical data: 

T v - 

~ constant ~ 12 keV (39) 

for light through heavy ground-state vector mesons (e q is the quark charge in units of e). 

This is the same as our data since the leptonic width is proportional to g 2 -^-. A final note 
on this result: a linear {Coulomb) potential implies a constant [linear) dependence on the 
meson mass for the decay constant, so our data is consistent with a superposition of the two 
potentials. 

In order to compare with experimental values of quarkonia we use results from Fig. [I] 
with masses nearest the physical J/ip and T masses. However, the masses are not tuned 
precisely to the experimental values (see Table |V|). For example, the (3 = 7.2 decay constants 
should be larger according to the derived afy cx y/M v a dependence. With this caveat, 
Figs. |2| and |3] show the comparison of our calculated decay constants with empirical values. 
Previous work on spin splittings would lead us to expect decay constants smaller than 
experimental values and this is what we find although the 10-15% underestimation (with 
the order v 2 relativistically corrected currents) may be a bit less than what one might have 
anticipated. As expected — since v 2 ~ .3 and v 2 ~ .1 — the shift from the order v 2 relativistic 
corrections in the currents is approximately three times as big for charmonium as compared 
to bottomonium. 

The main B c (spin averaged ground state ~ 6.35 GeV; see Table |TV|) results are shown in 
Fig. f|. The same information is in Table [VTJ. In Fig. |5| a comparison with previous lattice 



results [|T0| , |TT| for fs c combined with a partial list of other model results pl|-p5| is shown. 
Our final result, 420(13) MeV, is taken from the f3 = 7.4 data point (a ~ .16 fm). 

This is an initial study of decay constants made up of heavy quarks which systematically 
includes removal of the leading finite lattice spacing errors in the action, and measures 
the effects of the leading relativistic corrections in the vector and axial vector currents. The 
current matrix elements are matched with continuum QCD through one-loop in perturbation 
theory at lowest order in v 2 . After removal of the leading finite lattice spacing errors in the 
action we find these matching corrections to be small. This is shown in Tables [TI| and |TI: 
10% for the J/^; 5% for the T; and less than 1% for the B c . 
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APPENDIX A: GLUON PROPAGATOR AND 5V LA t FOR THE B c 



The gluon propagator follows from the quadratic part of Sq + Sgf- In the Feynman 
gauge (£ = 1) the gluon propagator is [19| (A is our gluon mass) 



D^ik) = D V(t {k) 



1 



1 



k 2 k 2 + X 2 



k\jk v + / ^{kgS ul / k v b \uj)k(jA. av {]z) 



(Al) 



with (ci 



•1/12) 



AUk)=A v Jk) = (l-S IJV )A(k)- 1 



(p) 2 - Cl P 2E^ 4 + fc 2 E K 



I V P / P T^p,,V P+P,V 



(A2) 



where 
A(k) 



k 2 - Cl Y,K 



fc 2 - Cl (P) 2 + E^ (^ 2 ) 3 + 2E^-^ 2 E^ 



The usual 



and 



fc M = 2sin(V2) 



(A3) 



(A4) 



P = 4^sin 2 (A; M /2) 
11=1 



(A5) 



definitions have been made. Also note that although we are working in the Feynman gauge 
here, is not diagonal in the Lorentz indices. 

The self-energy renormalization is performed as described by Morningstar in [I3J]. The 
vertex corrections follow in a similar manner. We will not write down all the Feynman rules 
here, but rather will do a particular example: the one-loop amputated vertex correction of 
the axial vector current at lowest order in v 2 in lattice NRQCD for the "free B c " system. 
In the notation of Eq. fl20|) , this is g 2 5Vi at . This is a standard vertex correction: the axial 
vector current creates a c and b which then exchange a gluon between them. The one- 
loop amputated vertex correction of the axial vector current at lowest order in v 2 in lattice 
NRQCD for the "free B c " system is 



4 f 7r 

g 2 5V lat =~g 2 



4 2 d 4 k D^{k)N t 



(A6) 



U (2vr)4 A c (-k)A b (k) ' 
the inverse propagator of the c is (F c and E c are defined below; n c is a stability parameter) 



12 



A c (-A;) = 1 - eM^)F^{k)E 2 c {k) ; (A7) 

the inverse propagator of the b is (F& and E\, are defined below; n& is a stability parameter) 

A b (k) = 1 - exp(-ih)F b 2nb (k)E^(k) ; (A8) 

D^ v is the above gluon propagator; and N^, a product of one gluon emission and absorption 
vertices, will be written below. There are no external momenta in Eq. (|A6|) because at lowest 
order in v 2 the three-momenta of the c and b are set to zero. Also recall that the mass has 
been subtracted from the energies as part of the definition of NRQCD. In the above inverse 
quark propagators 

F(k) = 1 - (A9) 
v ' Amn y ' 

and 

E (k) = 1 - Y -p- + 1*-L- , (A10) 
^ 48m 32nm 2 ' v ' 

for the particular c and b masses and stability parameters. The definitions 

fci = 2 sin(Aci/2) (All) 

and 

3 

k 2 = 4$]sin 2 (A; i /2) (A12) 
i=i 

have been made. The final piece is N^. It is defined as 

D^N^ = N T + N S + N ST , (A13) 

where the 'S' and 'T' subscripts signify 'spatial' and 'temporal' components respectively^ 
These three terms are 

(1) N T = D M (k)E b {k)E c {k)F b n "{k)F^(k) , 

(2) N s = £ D v {k)N v , 

and (3) N ST = £ D l4 (k) (N l4 + N 4i ) , (A14) 
i=i 

where 



Recall that our gluon propagator is not diagonal in Lorentz indices — hence the spatial-temporal 
terms. 
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Nij = sin(ifei/2) sin(%/2) 



1 



-£ c (fc)S< (0, fc) {1 + exp(^ 4 )F c ^(fc)} 



x 



2m c n c 

+2Kf(0, k) {l + esxp(ik 4 )Ff^(k)E t 
1 



0/ ) ^(^(0, k) {1 + exp(-i*4)l!?»(*)} 
+2^(0, k) {l + expHA^F^A;)^*;)} 



5 n (p,A;) = 5:F Q - 1 (p)F"- a (A;) 



a=l 



and 



/4 



1 



2m r n, 



sLn(ki/2)E c (k)S c nc (0, k) {1 + exp(^ 4 )F c ^(A;)} 



+2sin(ib i /2)iif?(0, k) H + exp(zfc 4 )F,f nc (A;)£, 



zexp(-ifc 4 /2)E fe (A;)F b n6 (fc) 



The specification of the integrand of Eq. ( \KBj ) is now complete. 
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APPENDIX B: ONE-LOOP AXIAL VECTOR CORRECTION 



This Appendix analyzes the origin of the '1/V terms in Eq. (p3|). First the NRQCD 
correction in Eq. ( |23a| ) is discussed, and then we move on to the QCD correction in Eq. ( |23b| ). 

The NRQCD correction in Eq. ( |23a|) : The integral of interest here is the complete one- 
loop amputated axial vector correction of continuum NRQCD. The integral is straightfor- 
ward and so we will not go into details here. Its integration is clearly described around 
Eqs. (14)-(16) of BF.Q Note that the result is entirely '1/V terms [compare with Eq. fllBp]. 

Now we discuss the correction in Eq. (|23b ) . The integral of interest is the infrared diver- 
gent piece of the one-loop amputated axial vector correction in continuum QCD (the other 
pieces are not discussed here because their integration is straightforward). The complete 
amputated vertex correction is written in Eq. (11) of BF.[] The term of interest is the first 
part of the first term in this equation: 



A 



IR 



dx 



l-x 



dyjj, 



2c 



d D k 



2 p ■ p' 



(2tt) d [k 2 - (xp 1 - yp) 2 + is] 



(Bl) 



The 4-momentum of the b is p 



p 2 + m 2 ,p), and in the center-of-mass frame the 4- 



momentum of the c is p' = {^jv 2 + m 2 , — p). After integrating over k, changing variables to 
s = x + y and t = x/s, and then integrating over s (as described in BF), the result is 



A 



IR 



2a, 
' 3tt 



■m b m c 



dt- 



At-ie 



1 



f< 2 



+ log. . 

e IR \A t -ie / 



l + 0(v 2 ) +0(e IR ), 



(B2) 



where the usual D = 4 — 2em replacement has been made; for the s-integral to converge, 
em has to be negative; jl 2 = 4iifi 2 /e' y ; and A t is given by 



A t = q 2 {t-t + ){t-t_ 



(B3) 



where 



and 



q 



(p + p) = (m b + m c ) + m b m c v + 0(v ) 
1 

q 2 



m h + p ■ p' ± J (p ■ p') 2 — mlm 2 



m b ± m red v 



m b + m c 

Recall l/m re ct = l/m b + l/m c . The origin of the terms is seen by rewriting 



(B4) 



(B5) 



m b m c 



m b m c 



A t - is q 2 {t + - t-) \t - t H 



IE 



t — t_ + ie 



(B6) 



7 In this Appendix we will take "Ref. |§" — > "BF". 

8 BF's notation is used in this Appendix. To convert to our notation take A — ► g 2 5V and 



V(4vr). 
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and noting that 



m b m c 



t. 



2v 



+ 0{v) 



Note that the 0(v 2 ) terms of Eq. (|B5| ) cancel in the difference: 



U 



t. 



2m red v 
m b + m c 



+ 0{v 3 ) . 



To proceed with the integration, we find Lewin's book on polylogarithms 
More specifically, we use the following three equations: 



u \og( a + bt) , 1, f ae - bc\ , fc + eu\ 1 
— — dt = - log log - - Li 2 

oc + er e \ e / V c / e 



6(c + eit) 



1 T . / 5c ' 
+ - Li 2 

e \bc — ae 



be — ae 
be — ae ^ . 



(B7) 

(B8) 
helpful. 



(B9a) 



and 



Li 2 (x ± ie) x = ^- - ^ log 2 (x) ± Z7r log(x) 



7T 1 

Li 2 (x) ^ -— - - log 2 (-x) 



\x) (x) (z) 



2 2 



+ 



+ 



3 2 



+ 



(B9b) 



(B9c) 



where a, 6, c, and e may be complex, but x is real. 142(2) is the dilogarithm function nicely 
elucidated in Ref. |2T| . 

The integration of Eq. ( |B2| ) is now straightforward with result (as v — > 0) 



A 



7\R 



2a, 
3vr 



— -2 + 21og^^- 

e IR m b + m c 



m red m red 
2 log 



7T 27T 



2 log 



2m red v \ 



2 log 



(BIO) 



where recall log(/i 2 ) = log(yU 2 ) + log(47r) —7. Note that all the l l/v' terms of the full answer 
[Eq. ( |23b|) 1 are here, as was to be shown. The imaginary 'l/e/i?' term is also here, but there 
appears to be an extra real 'l/e^' term. The fact that this is not an "extra term" is seen 
after including the wave function renormalization [Eq. (10) of BF]: 



2a, 
3tt 



46 



uv 



1 3 m Q 
h - log — 

2e IR 2 fj, 



(Bll) 



both the b and c contribute one of these factors which results in the cancelation of the 
real 'l/e^' term of Eq. ( [B10|) . After including both of these wave function renormalization 
factors and the complete amputated vertex correction, the result is Eq. ( |23b| ) as mentioned 
in the body of the paper. 
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FIG. 1. Meson mass dependence of the order v 2 relativistically corrected vector decay con- 
stants of quarkonia including the one-loop perturbative matching. The dotted line is proportional 
to V / W a. 
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FIG. 2. The calculated vector decay constant of charmonium (including the perturbative 
matching) at different values of the lattice spacing with and without the relativistic corrections in 
the current. The dotted line is the experimental result. 
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FIG. 3. The calculated vector decay constant of bottomonium (including the perturbative 
matching) at different values of the lattice spacing with and without the relativistic corrections in 
the current. The dotted line is the experimental result. 
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FIG. 4. Scaling behavior of the B c decay constant (including the perturbative matching) with 
and without the relativistic corrections in the current. 
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FIG. 5. A comparison with previous work. The filled square is this work: 0(1, f 2 ); filled 
triangle is this work: 0(1); and the open circles are other models labeled by their respective 
reference. 
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TABLES 



TABLE I. Comparison of perturbative (one loop) and nonperturbative Landau links using the 
boosted definition of the coupling [see Eq. (|lC|)], Recall uq = 1 — chUq . 



p 


u = (\ReTrU^) 


a = 10/(471-^4) 


(2) 

nonpert 


pert 


7.2 


.805 


.2632 


.7409 


.7503 


7.3 


.817 


.2447 


.7480 


.7503 


7.4 


.8286 


.2281 


.7513 


.7503 



TABLE II. One-loop perturbative lattice NRQCD matching factors of quarkonia. Z matc h is 
defined in and above Eq. (|2"§|). The results of the middle column are to be multiplied by 4/[3(27r) 4 ]. 
The statistical error estimate from VEGAS is less than one percent. 



ma[n] 


SVlat + SZg 


•^match 


1.4[3] 


-43.9335 


0.9141 


1.5[3] 


-42.3631 


0.9037 


1.6[3] 


-45.3303 


0.9049 


2.1[2] 


-49.3386 


0.9274 


2.2[2] 


-49.8716 


0.9235 


2.3[2] 


-50.3880 


0.9192 


2.8[2] 


-52.8248 


0.9359 


2.9[2] 


-53.2263 


0.9323 


3.0[2] 


-53.7950 


0.9288 


3.5[2] 


-57.0103 


0.9462 


3.6[2] 


-57.5486 


0.9437 


3.7[2] 


-58.3099 


0.9416 


4.2[2] 


-61.5731 


0.9574 


4.3[2] 


-62.3861 


0.9564 


4.4[2] 


-63.3444 


0.9558 



TABLE III. One- loop perturbative lattice NRQCD matching factors of the B c . ^ ma tch is defined 
in and above Eq. (p8[). The results of the middle column are to be multiplied by 4/[3(27r) 4 ]. The 
statistical error estimate from VEGAS is less than one percent. 



m c a[n], mi,a[n] 


+ + 52^/2 


•^match 


1.4[3],3.8[2] 


-47.5372 


1.0048 


1.5[3],4.4[2] 


-47.2181 


1.0096 


1.6[3],5.0[2] 


-49.3620 


1.0213 
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TABLE IV. Mp Crt versus M k i n for runs tuned to 6.35 GeV for the spin averaged B c . The 
statistical error estimate of M pert from VEGAS is less than one percent. 






a(fm) 


m c a[n], niba[n] 


M kin a 


Mp ert a 
w/o tad imp w 


tad imp 


M kin (GeV) 


7.2 


.204(1) 


1.6[3], 5.0[2] 


6.60(21) 


6.91 


6.92 


6.37(21) 


7.3 


.185(2) 


1.5[3], 4.4[2] 


5.94(16) 


6.15 


6.28 


6.32(18) 


7.4 


.163(3) 


1.4[3], 3.8[2] 


5.26(18) 


5.42 


5.67 


6.37(24) 



TABLE V. Mp Cr t versus M k ; n for runs nearest the charm and bottom regions respectively. The 
statistical error estimate of Mp er t from VEGAS is less than one percent. 



Mpert a 



p 


a(fm) 


ma[n] 


M kin a 


w/o tad imp 


w tad imp 


M kin (GeV) 


7.2 


.240(3) 


1.6[3] 


3.37(6) 


3.17 


3.89 


2.77(6) 


7.3 


.205(3) 


1.5[3] 


3.20(5) 


2.96 


3.69 


3.08(6) 


7.4 


.178(3) 


1.4[3] 


3.00(6) 


2.80 


3.55 


3.32(9) 


7.2 


.198(2) 


4.4[2] 


9.07(21) 


9.27 


8.78 


9.05(23) 


7.3 


.174(3) 


4.3 [2] 


8.74(20) 


9.06 


8.63 


9.91(28) 


7.4 


.151(2) 


3.5[2] 


7.05(20) 


7.35 


7.19 


9.20(29) 



TABLE VI. Vector decay constants of quarkonia (after including perturbative matching). 



fv (MeV) 



p 


ma[n], M kin a 


a(fm) 


w/o v 2 rel cor 


w v 2 rel cor 


7.4 


1.4[3], 3.00(6) 


.178(3) 


411(8) 


364(9) 


7.3 


1.5[3], 3.20(5) 


.205(3) 


396(7) 


352(8) 


7.2 


1.6[3], 3.37(6) 


.240(3) 


375(6) 


335(8) 


7.4 


2.1[2], 4.33(10) 


.164(2) 


492(8) 


464(10) 


7.3 


2.2[2], 4.50(7) 


.188(2) 


474(6) 


446(8) 


7.2 


2.3[2], 4.69(8) 


.220(2) 


440(6) 


415(8) 


7.4 


2.8[2], 5.66(15) 


.156(2) 


561(10) 


541(12) 


7.3 


2.9[2], 5.84(11) 


.179(2) 


533(8) 


515(10) 


7.2 


3.0[2], 6.06(13) 


.208(5) 


495(13) 


477(14) 


7.4 


3.5[2], 7.05(20) 


.151(2) 


619(12) 


605(14) 


7.3 


3.6[2], 7.26(15) 


.174(2) 


581(9) 


567(11) 


7.2 


3.7[2], 7.52(17) 


.201(2) 


538(8) 


524(10) 


7.4 


4.2[2], 8.51(29) 


.149(2) 


662(14) 


651(15) 


7.3 


4.3[2], 8.74(20) 


.174(3) 


608(12) 


597(13) 


7.2 


4.4[2], 9.07(21) 


.198(2) 


567(9) 


557(11) 
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TABLE VII. B c decay constants (after including perturbative matching). 






a(fm) 


f Bc (MeV) 

w/o v 2 rel cor 


w v 2 rel cor 


7.2 


.204(1) 


462(8) 


407(11) 


7.3 


.185(2) 


463(8) 


406(10) 


7.4 


.163(3) 


478(12) 


420(13) 
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